Projective Quantum Monte Carlo Method for the Anderson Impurity Model 
and its Application to Dynamical Mean Field Theory 
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We develop a projective quantum Monte Carlo algorithm of the Hirsch-Fye type for obtaining 
ground state properties of the Anderson impurity model. This method is employed to solve the 
self-consistency equations of dynamical mean field theory. It is shown that the approach converges 
rapidly to the ground state so that reliable zero-temperature results are obtained. As a first ap- 
plication, we study the Mott-Hubbard metal-insulator transition of the one-band Hubbard model, 
reconfirming the numerical renormalization group results. 
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In recent years, we have seen a revival of interest in 
Kondo-like physics, in particular in quantum dot systems 
0, adatoms on surfaces 0, and in connection with the 
dynamical mean field theory (DMFT) 0, Q . The under- 
lying microscopic model, the Anderson impurity model 
(AIM), can be solved exactly by Bethe ansatz |5| for the 
special case of a linear dispersion relation. But in general, 
such a solution is only possible numerically. Two sophis- 
ticated numerical techniques have been established in the 
literature: the numerical renormalization group (NRG) 

3 and the Hirsch-Fye quantum Monte Carlo (HF-QMC) 

JJ method. 

In his pioneering work 6], Wilson applied the renor- 
malization group concept to the Kondo problem: High 
energy degrees of freedom are systematically integrated 
out and one obtains an excellent description of the low 
energy degrees of freedom, which was not possible before 
since perturbation theory for the AIM breaks down at 
low temperatures and energies. Confirming earlier scal- 
ing ideas of Anderson these NRG calculations show 
the formation of a low-energy Kondo singlet between the 
AIM localized and conduction electrons and an associate 
Abrikosov-Suhl resonance in the spectrum. 

The same physics can be described by the HF-QMC 
algorithm. This method solves the AIM numerically by 
discretizing the imaginary time interval [0, j3\ (f3 = 1/T; 
inverse temperature) into steps of size At and by map- 
ping the interacting AIM onto a sum of non-interacting 
models via the Hubbard-Stratonovich transformation. 

The two approaches require the extrapolation of a dis- 
cretization parameter, the logarithmic energy mesh pa- 
rameter A -> 1 (NRG) and At -> (HF-QMC), and have 
some limitations: The NRG effort grows exponentially 
with the number of localized AIM orbitals M, restrict- 
ing the algorithm to M < 2, which prevents the usage 
of NRG for more realistic calculations where one would 
like to include, e.g., M — 5 or 7 orbitals for 3d and 4/ 
systems, respectively. The HF-QMC on the other hand 
scales like M 2 (f3/Ar) 3 and, hence, quickly becomes too 
expensive in CPU time if temperature is decreased. This 
is a severe restriction since interesting many-body physics 
often occurs at low temperatures. 



In this Letter, we introduce a new projective quan- 
tum Monte Carlo (PQMC) method for studying the zero- 
temperature AIM in the thermodynamic limit. Similar 
to projective approaches in lattice models |9j, it relies 
on the idea that ground state properties are more effi- 
ciently obtained by filtering out - via projection along 
the imaginary time axis - the ground state wave func- 
tion from a suitably chosen trial wave function. For the 
AIM, however, an efficient algorithm needs to be of the 
Hirsch-Fye type, in contrast to the lattice case 9J. In a 
second step, we incorporate this algorithm as an impu- 
rity solver within DMFT. In passing, we provide for a 
mandatory confirmation of the T=0 DMFT scenario for 
the Mott-Hubbard transition. 

PQMC algorithm - Our starting point is the AIM 
Hamiltonian 
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Here, create (annihilate) impurity electrons with 

spin a which interact via U and have a level energy e/; 
v,f a — /ct/ct- These orbitals hybridize via V k with a 
conduction band (ct), having a dispersion e, . 

Let us consider a trial wave function |\&t) which is 
non-orthogonal to the (non-degenerate) ground state 
|\&o) of the AIM. By means of the 8- projector ~ 
cxp(— 0Haim/2), we can filter out \^q) from |^t) and 
calculate ground state expectation values of an arbitrary 
operator O: 



(O)o = 
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We now relate the projection formula to an artifi- 
cial finite temperature (1//3) problem, using an auxiliary 
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Hamiltonian Ht whose ground state is \^>t)' 
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With this trick, the r.h.s. of Eq. Q is in a suitable form to 
apply (or rederive) the Hirsch-Fye algorithm at finite val- 
ues of (3. In the following, we will point out the necessary 
steps in short, more details will be published elsewhere 

u. 

We discretize imaginary time into Trotter slices of 
length At and decouple the Coulomb interaction via the 
Hubbard-Stratonovich transformation. Let us now re- 
strict to single-particle trial wave functions \^>t) and 
Hamiltonians Ht- Then, the interaction (and the 
Hubbard-Stratonovich field) is zero for the (3 part, 
and we can combine the propagation of this part 
(exp[— (3Ht\) with the H part of the first 9 Trotter 
slice (exp[— AtHq])- The limit (3 — > oo can now be 
taken analytically, leaving a problem on the interval 
[0,9] discretized into L = 9/ At steps. As the HF- 
QMC, our PQMC algorithm deals with L x L matri- 
ces for the f-electron Green function on this interval. 
The updating equations for these Green functions (if a 
Hubbard-Stratonovich field is changed) are the same as 
that of HF-QMC. Choosing Ht = H a , also the start- 
ing Green function (with all Hubbard-Stratonovich fields 
zero) has a similar structure: In HF-QMC it is the 
finite-temperature non-interacting Green function. In 
our PQMC algorithm it is - because of the extra (3 — > oo 
projection - the zero-temperature non-interacting Green 
function Qq(t,t') truncated to < t,t' < 9. But with 
this different object (initial Green function), the same 
program code can actually be used. 

Our PQMC algorithm allows the measurement of the 
Green function G (r, r') but also of arbitrary correlation 
functions, (^-projected from the slater determinant of the 
non-interacting AIM as a trial wave function \^t)- For 
a meaningful PQMC calculation, we should measure all 
correlation functions only on the central C time slices. 
Then V = (L — C)/2 time slices on the right and left 
side of the measuring interval serve for projection. In 
the following, QMC results are understood in this sense. 

DMFT self- consistency implementation - We now im- 
plement this algorithm as an impurity solver in the 
DMFT self-consistency cycle. Let us consider the half- 
filled Hubbard model 
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where denotes the restriction to nearest neighbor 

hopping in d dimensions. We use the Bethe lattice with 
semicircular density states N(E) — -^w?\/W 2 /4: — E 2 . 
In the following, units are normalized to a bandwidth 
W = 4. 



DMFT maps the Hubbard model onto an AIM with 
non-interacting Green function 



Go{iw n ) = [G 1 (iuj n ) + £(iw)]~ 



(6) 



where u n denote Matsubara frequencies. The Green 
function G(ioj n ) of the AIM has to be determined self- 
consistently together with the self energy T,(ioj). Both 
are connected by the k-integrated Dyson equation 



G(iuj n ) = [dE- 
J * 



N(E) 



(7) 



iu — E — Y,{iu) n ) 
Thereby, the DMFT self-consistency cycle is as follows: 

l^yiiOn) ~ G(lLU n ) =5 y (lLU n ) — > G{T) =3 h{lLU n ) ■ ■ ■ 

For the difficult task, i.e., calculating the zero tempera- 
ture Green function G(r) of the AIM for a given Go(t) 
with r £ [— 9, 01, we use the PQMC algorithm, instead 
of HF-QMC Since the PQMC Green function G(t) 
is only defined for r < CAt, we first extrapolate G (t) 
to large r's before Fourier transforming to Matsubara 
frequencies in the above self-consistency cycle. To this 
end, we employ the maximum entropy method, yield- 
ing the spectral function A (oj) which allows to calculate 
G(iLu n ) = J duj gij~±; at any frequency iu n . It is impor- 
tant to note that, in the end, this procedure only serves 
as a fit in imaginary time. Hence, we do not suffer from 
the usual problem - determining A(uj) from G(r) is ill- 
conditioned. Only for very small iw n , corresponding to 
the extrapolation in imaginary time, this procedure even- 
tually breaks down. 

Results - We will now use our PQMC method to study 
the interaction driven Mott-Hubbard metal-insulator 
transition that occurs at half filling. First, in Fig. we 
illustrate the power of the PQMC method in direct com- 
parison with finite-T HF-QMC results. We choose equal 
values of inverse temperature f3 and total projection time 
9, having an identical number of Trotter time slices for 
both calculations. Then, our results for (3 = 9 are ob- 
tained with a comparable numerical effort. For U = 4.8 
and at low T, the DMFT equations have two solutions, 
where we follow the metallic solution (as soon as it ex- 
ists). In the case of DMFT calculations with finite-T 
HF-QMC, the double occupancy per site D = {n^ni) of 
Fig- B starts at high temperatures ((3 = 10) in the insu- 
lating phase and only after the transition to the metallic 
phase around f3 = 40 convergence to the ground state 
value sets in. In contrast, the PQMC double occupancy 
convergences much faster and extrapolates almost lin- 
early to the T=0 value. For example, the 9 = 10 PQMC 
result is even closer to the T=0 result than the HF-QMC 
double occupation at (3 — 30. Since the numerical effort 
grows cubically with 9 or (3 this means that the PQMC 
is roughly 27-times faster at the same accuracy. We at- 
tribute this behavior to the absence of thermal fluctua- 
tions in PQMC. 
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FIG. 1: DMFT double occupancy D — (n^n{) as a func- 
tion of (i) temperature T = 1//3 (conventional HF-QMC; cir- 
cles) and (ii) the inverse projection parameter 1/8 (PQMC, 
squares). For 8 = f3 the numerical effort is comparable, but 
the PQMC results converge much faster and almost linearly 
to the T — value. Here and in the following, the bandwidth 
W = 4 sets the energy unit. 



In the following, we present our results for the Mott- 
Hubbard metal-insulator transition, starting with the 
double occupation in Fig. [21 All calculations are done 
at a fixed At = 0.2. Upon increasing U, we systemati- 
cally follow the metallic solution until at U C2 w 6.0 the 
metallic quasiparticle peak disappears. For the insulating 
solution on the other hand, we decrease the interaction 
until the insulator becomes unstable below U Cl « 5.0. 
Within the uncertainty of roughly 0.1-0.2, our coexis- 
tence region U C1 < U < U C2 compares well with NRG 
results (U C1 = 4.78, U C2 = 5.88 [H EH). Already the 
(3 = 20 results (circles) for the double occupancy yield 
a good estimate for the low temperature coexistence re- 
gion. In order to resolve the increasingly sharp quasipar- 
ticle peak in the vicinity of the critical point we also used 
longer projection times 8 and the thick line represents a 
linear extrapolation for three values 8 = 20,30,40. Also 
plotted is the result from a 10th order perturbation ex- 
pansion H3 in t/U, which agrees well with our results 
for the insulator. Our results are consistent with a linear 
dependence (D mct — D ins ) oc (U — U C2 ) which may be inte- 
grated to a groundstate energy (E mct —E ins ) oc (U—IL 2 ) 2 
corresponding to a second order transition at U C2 16]. 

Fig. |31 plots the local spin susceptibility Tx\ oc (0) = 
lim \ Jq (S z (t) S z ) dr, and actually represents the first 

DMFT analysis of this quantity through the zero temper- 
ature Mott-Hubbard metal-insulator transition (at finitc- 
T it was studied in EH)- In the top panel, we observe the 
local moment behavior of the insulating solution without 
any sign of criticality at U C2 . On the other hand ap- 
proaching U C2 from the metallic side, xi oc is expected to 
diverge, reflecting the critical behavior of a Fermi liquid 
with a vanishing quasiparticle weight Z 16]. Therefore 
the second order character of the phase transition at U C2 
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FIG. 2: Double occupancy D — (n-j-nx) as a function of inter- 
action U for both metallic and insulating DMFT solutions. In 
the interval 5.0 « U C1 < U < U C2 ~ 6.0 both solutions coexist 
(shaded region). Circles: 8 — 20; thick line: extrapolation to 
8 — » oo; dashed line: strong coupling expansion of (l4|. 
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FIG. 3: Static local spin susceptibility Txioc(O) vs. U for the 
insulating (upper panel) and metallic DMFT solution (lower 
panel). For the metal, Xioc(0) diverges as we approach U C2 
from below. Since we plot the truncated integral xSci the 
divergence is cut off at Xioc = 37.4 (star). 



can only be observed for U < U C2 . For the metallic phase, 
we integrate the spin correlation function up to a max- 
imal cutoff C: x? oc (0) = fa (S z (r) S z ) dr. The thick 
line is an extrapolation of X\ oc f° r a nxe d cutoff C = 40 
to 8 — > oo, i.e., to the groundstate value. At the critical 
point, we see indications for a divergence of xi oc which 
for the cutoff quantity X\o C means approaching the star 



(Xioc = 37.4, the value of the insulating solution). 

Finally, we present in Fig. 0] the spectra at U = 5.9, 
close to the phase transition. The insulating solution 
(8 = 20) has a pronounced charge gap of around 1.2 be- 
tween two Hubbard bands. The metallic solution shows 
an additional quasiparticle peak which is pinned to the 
non-interacting value, a consequence of Fermi liquid the- 
ory at T = 17]. The quasiparticle peak for a fixed 
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FIG. 4: Spectral functions A (ui) for the insulating (line) and 
metallic solution (dashed) close to the phase transition at U = 
5.9). Inset: Quasiparticle weight Z vs. U for 6 — 20 (circles) 
and extrapolated to 9 — > oo (thick line) , compared with NRG 
E3 (dashed line). 

9 = 40 has width of 0.05 (Note, that the low energy 
accuracy is limited by 1/9 and the statistical quality 
of QMC data.). We determine the quasiparticle weight 

by Z = (l — — ^ UJl - J from the first Matsubara fre- 
quency. In principle, c_>i can be chosen to be arbitrarily 
small, but controlled values for E (ioj n ) are restricted by 
the projection length 1/9 and the maximum entropy ex- 
trapolation fit. The inset in Fig. 0] shows that at 9 = 20 
the weight Z is still too large. But after extrapolating Z 
from 9 = 20,30,40 to 9 — ► oo, we find agreement with 
the NRG results [12. 

Summary - We have introduced a novel projective 
(PQMC) algorithm for impurity models, which is on the 
technical front closely related to the standard Hirsch-Fye 
algorithm. The numerical effort for both algorithms is 
identical, but the convergence to the ground state is dra- 
matically improved by PQMC. Using PQMC as an im- 
purity solver for DMFT, we were able to obtain accurate 
zero temperature results for the Mott-Hubbard metal- 
insulator transition. There has been a strong controversy 
about the nature of this Mott-Hubbard transition within 
DMFT, in particular, whether two solutions (metallic 
and insulating) coexist in an interval U c \ < U < U C 2 
or not. While the coexistence is by now generally ac- 
cepted, NRG calculations represent the only unas- 
sailable zero-temperature confirmation of the coexistence 
scenario Since other numerical results at T = 0] 
and analytical arguments |19| contradict this scenario, an 
independent numerical study is certainly in order. Using 
PQMC, we have provided for such a mandatory reconfir- 
mation, and beyond it we presented results for the double 
occupancy and local susceptibility. 

Outlook - In recent years, there have been two main 
lines of DMFT developments: towards realistic calcula- 
tions such as LDA+DMFT and towards the inclu- 



sion of more and more (non-lo cal) correlations by means 
of cluster DMFT calculations |2lJ . These approaches re- 
quire additional orbitals and DMFT sites, respectively, 
making NRG prohibitively CPU-expensive since the ef- 
fort scales exponentially with more orbitals and sites. 
In contrast, our PQMC algorithm only scales quadrat- 
ically or cubically, and hence opens the door to accessing 
low temperatures within realistic approaches or model- 
oriented calculations with non-local correlations. 
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